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ABSTRACT 

Based on our high resolution, two-dimensional hydro dynamical simulations, we pro- 
pose that large cavities may be formed by the nonlinear development of the combined 
thermal and gravitational instabilities, without need for stellar energy injection in a 
galaxy modeling the Large Magellanic Clouds (LMC). Our numerical model of the star 
formation allows us to follow the evolution of the blastwaves due to supernovae in the 
inhomogeous, multi-phase, and turbulent-like media self-consistently. Formation of kpc- 
scale inhomogeneity, such as cavities, observed H I map of the LMC, is suppressed by 
frequent supernovae (average supernova rate for the whole disk is ~ 0.001 yr _1 ). However 
the supernova explosions are necessary for the hot component (T g > 10 6-7 K). Position- 
velocity maps show that kpc-scale shells/arcs formed through the nonlinear evolution in 
a model without stellar energy feedback has similar kinematics to explosional phenom- 
ena, such as supernovae. We also find that dense clumps and filamentary structure are 
formed due to a natural consequence of the non-linear evolution of the multi-phage ISM. 
Although the ISM in a small scale looks turbulent-like and transient, the global struc- 
ture of the ISM is quasi-stable. In the quasi-stable phase, the volume filling factor of the 
hot, warm, cold components are ~ 0.2, ~ 0.6, and ~ 0.2, respectively. We compare the 
observations of H I and molecular gas of the LMC with the numerically obtained H I and 
CO brightness temperature distribution. The morphology and statistical properties of 
the numerical H I and CO maps are discussed. We find that the cloud mass spectrum of 
our models represent a power-law shape, but their slopes change between models with 
and without the stellar energy injection, and also the slope depends on the threshold 
brightness temperature of CO. 

Subject headings: ISM: structure, kinematics and dynamics — galaxies: structure - 
individual: LMC — method: numerical 
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1. INTRODUCTION 

The topology of the neutral interstellar medium 
(ISM) can be studied in great detail by the spatial 
and velocity structures in the neutral H I gas. A 
recent high-resolution H I survey of the Large Magel- 
lanic Cloud (LMC) reveals that the structure of the 
neutral atomic inter stellar gas is dominated by nu- 
merous holes a nd shells as well as complex filamen- 
tary structure (Kim et al. 1998). These features are 



commonly seen in recent high-resolution H I images of 
nearby galaxies obtained with radio synthesis intcrfer- 



ometers (Deul & den Hartog 199C; 


Puchc et al. 1992: 


staveley-Smith et al. 1997|; 


Walter & Brinks 199£; 


stanimirovic et al. 1999). In general, the shell-like 



and hole structure seen in H I has been understood 
as the cumulative effect of stellar winds from massive 
stars and supernova explosions evacuating the cool 



ISM ( 


Tenorio-Tagle 1988 




van dcr Hulst 1996 


; Oey 


199G 




Oey & Clarke 1997 


). However, the extensive 



study of H I shells in the LMC shows that there is a 
relatively weak correlation between the H I shells and 
the ionized gas traced out by the H II regions and 
H II filaments (Kim et al. 1999). Furthermore, the 



correlation between the H I shells and 122 OB stellar 



assoc iations in the LMC is not very tight ( Kim & Chu 
2000 ). Rhode et al. (1999) claimed that there is no 
remnant star clusters at the center of the H I holes in 
Holmberg II, and it is inconsistent with the SNe hy- 
pothesis. Moreover an energy source generating the 
kpc-scale supergiant H I holes is a puzzle. These is- 
sues raise an interesting question about whether H I 
shells/holes have been formed by the interaction be- 
tween stars and the ISM or not. 

Recent hydrodynamical simulations by Wada & 
Norman (1999) demonstrate that a gravitationally 
and thermally unstable disk, which models an ISM 
in galaxies, can generate the cold, dense clumps and 
filaments surrounded by hot, diffuse medium. They 
show that porous structure is a natural consequence 
of the non-linear evolution of the ISM. This result 
strongly suggests that some fraction of H I shells, su- 
pershells or holes seen in galaxies does not relate to 
the interaction between stellar activities and the ISM. 

The other important component of the ISM is the 
molecular gas. Molecular clouds are potential sites of 
star formation, but their formation mechanism, and 
the relationship between their evolution and star for- 
mation have been poorly understood. Dense molec- 
ular hydrogen is traced by the rotational transition 



of CO, and a recent high resolution survey of 12 CO 
(J=l-0) in the LMC with NANTEN, which is 4m 
millimctcr-wavc telescope at Las Campanas Obser- 
vatory, established a comprehensive view of the giant 
molecular clouds in the LMC ( Fukui et al. 1999| ) . 
Since the gas clumps, whose size is typically 10 — 100 
pc, in the simulations by Wada & Norman (1999) are 
dense (n > 1000 cur 3 ) and cold (T = 10 - 100 K), 
they are counterparts of the observed giant molecular 
clouds in the LMC or other galaxies. 

Using numerical simulations, Feitzinger et al. (1981) 
and Gardiner, Turfus, & Putman (1998) have studied 
global dynamics of the ISM and star formation in the 
LMC. Although the numerical methods they used are 
different, the stochastic self-propagating star forma- 
tion model (e.g. Seiden & Gerola 1984) in Feitzinger 
et al. (1981) and sticky particle method in Gardiner, 
Turfus, & Putman (1998), are both phenomenolog- 
ical concerning the structure and dynamics of the 
ISM and star formation processes. Unfortunately the 
spatial resolution (100 pc in Feitzinger et al. 1981) 
and mass resolution (6xl0 4 M Q , which corresponds 
to ~ 200 pc for n ~ 1 cm -3 and H ~ 100 pc, in Gar- 
diner, Turfus, & Putman 1998) in these simulations 
are not good enough for comparison between models 
and the recent high resolution observations (~ 15 pc 
for HI (|Kim et al. 1998|) and - 40 pc for CO (J=l-0) 



( |Fukui et al. 1999j )). 

In this paper, we apply the numerical scheme 
used in Wada & Norman (1999) to a LMC-type 
model galaxy, and we conduct two-dimensional hy- 
drodynamical simulations of the multi-phase ISM in 
a LMC-like model galaxy, taking into account self- 
gravity of the gas, radiative cooling and various heat- 
ing processes, such as supernova explosions. High 
spatial resolution (7.8 pc) and a modern hydrody- 
namical scheme allow us to model the star formation 
and its feedback less phenomenologically, and there- 
fore it needs less assumptions than the previous semi- 
analytic and numeric al approaches (e.g. a review by 
Shore fc Ferrini 1995 ). In contrast to the model used 
in Wada & Norman (1999), the rotation curve as- 
sumed here is nearly rigid as suggested in the LMC 
or other LMC-type dwarf galaxies. We derive the H I 
and CO brightness map from the simulations. Then 
we compare these simulation results with the H I and 
CO observations of the LMC. 

In §2, we describe our numerical method and mod- 
els. In §3, the numerical results with and without star 
formation are discussed on morphology and statisti- 
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cal structure of the ISM, then they are compared with 
the observations. Position- Velocity diagrams derived 
from the numerical results are also discussed. Com- 
parison with past numerical simulations and other im- 
plications are discussed in §4, and conclusions are pre- 
sented in §5. 

2. NUMERICAL METHOD AND MODELS 

Taking into account the multi-phase and inho- 
mogenous nature of the ISM is crucial for realistic 
simulations of global star formation in galaxies. Stars 
are formed preferentially in cold molecular gas, and 
supernovae produce low density, high temperature 
gas. Interaction between such different phases are 



also impornant (e.g. 


McKee & Ostrikcr 1977; 


Ikeuchi, 


Habe, & Tanaka 1984 


) . In order to simulate dynamics 



of the inhomogeous, multi-phase interstellar matter 
and star formation, our numerical method has fol- 
lowing features. (1) High spatial resolution (7.8 pc), 
for diffuse gas to high density gas. We use an Eu- 
ler mesh code with 1024 2 Cartesian grid points. (2) 
The simulations are global in order to study effects 
of galactic rotation and other kpc scale phenomena 
to the local structure of the ISM. (3) Self-gravity of 
the gas is calculated. (4) Radiative cooling for gas 
whose temperature is between 10 s K and 10 K is 
taken into account. (5) Various heating processes are 
included. Here we assume photoelectric heating due 
to dust grains and UV radiation, SN explosions, and 
stellar wind from massive stars. (6) High numerical 
accuracy for shocks is achieved with a modern hydro- 
dynamical scheme. This is crucial, because the ISM 
is usually supersonic, and SNe produce strong shocks. 

We solve the following hydrodynamical equations 
and the Poisson equation numerically in two dimen- 
sions to simulate the evolution of a rotating gas disk 
in a stellar potential. 



dp 
at 



v • (H = o, 



(i) 



dv , Vp ^ x 

— + (v ■ V » + — + V$ cxt 
ot p 



V<£ sg = 0, (2) 



dE 1 

— + - V • [{ P E + p)v] = r uv + - pA(T g ), (3) 
ot p 



where, p,p,v are the density, pressure, and ve- 
locity of the gas, and the specific total energy E = 
\v\ 2 /2 + p/("f — l)p, with 7 = 1.4. We assume a 
time-independent external potential <I> cx t oc v 2 /(R 2 + 
a 2 ) 1//2 , where a = 2.5 kpc is a core radius of the po- 
tential and v c = 63 km s" 1 is the maximum rotational 
velocity, and they are determined to mimic the rota- 
tion curve of the LMC ( |Kim ct al. 1998| ). Since the 
total gas mass is about 10 % of the total dynami- 
cal mass (see below), the effect of self-gravity of the 
gas to the rotation curve is not significant. In fact 
the rotation curve after the system evolves is close to 
the rigid rotation (see §3.3 and the position-velocity 
diagrams in Fig. 10). 

We also as sume a cooling function A(T g ) (10 < 



T g < 10 8 K) (|Spaans fc Norman 19970 . The cooling 



processes taken into account are, (1) recombination 
of H, He, C, O, N, Si and Fe, (2) collisional excita- 
tion of HI, CI- IV and OI- IV, (3) hydrogen and helium 
bremsstrahlung, (4) vibrational and rotational excita- 
tion of H2 and (5) atomic and molecular cooling due 
to fine-strucure emission of C, C + and O, and ro- 
tational line emission of CO and H2. As a heating 
source, Tuvi we assume a uniform UV radiation field 
(Gerritsen & Ickc 1997), which is normalized to the 



V 2 $ sg = iirGp, 



(4) 



local interstellar value, and photoelectric heating by 
grains and PAHs. The bulk of the heating for the 
10 2 — 10 4 K gas is provided by photo-emission of UV 
irradiated dust grains (Bakes & Ticlcns 1994). 

We consider two feedback effects of massive stars 
on the gas dynamics, namely stellar winds and super- 
nova explosions, although our results show that the 
former is less effective than the latter. We first iden- 
tify cells which satisfy criteria for star formation. The 
criteria are a surface density threshold (£ g )ij > S c 
and a critical temperature {T g )i.j < T c below which 
star formation is allowed. The surface density is de- 
fined as £ g = 2Hp, where H is the scale height and 
assumed to be constant (100 pc). In these simulations 
we take S c = 40M© pc" 2 and T c = 15 K (model 'Star 
Formation 1', hereafter 'SF1') and ten times larger 
threshold density, E c = 400M Q pc" 2 and T c = 15 
K (model SF2). These criteria are chosen to sat- 
isfy the condition, Lj < 0.1 A or Lj < 0.01 A, where 
Lj and A are the Jeans length and the size of each 
cell. The volume filling factors of the star-forming 
cells to the total volume are typically ~ 5 x 10 -3 
and ~ 5 x 10 -4 for model SF1 and SF2, respectively. 
Namely the star forming sites are cold, clumpy re- 
gions. We also assume the star-forming criteria must 
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last during 10 5 yr for each cell before the star for- 
mation is initiated. Since the spatial resolution is fine 
enough compared to the GMC size, we do not have to 
assume any global criteria for gravitational instability 
to identify star forming sites. Assuming the Salpeter 
IMF with m u = 120M Q and mi = 0.2M Q , we cre- 
ate test particles representing massive stars (> 8M Q ) 
in the star forming cell. The kinematics of the test 
particles in the external potential and the self-gravity 
potential of the gas are traced by the second-order 
time-integration method. The stars (test particles) 
inject energy due to stellar winds d uring their lifetime 
wh ich is approx imately ~ 10 7 yr ( Lcithcrcr, Robert, 
& frisson 1992 ). When one of stars represented by 
a test particle explodes as supernova, an energy of 
10 51 ergs is injected into the cell as thermal energy 
where the test particle is located at that moment. 
The cooling procedure is not used for such cells, but 
the cells adjoining the supernova cell are treated nor- 
mally. Evolution of the supernova remnant is very 
dependent on its environment. In contrast to past 
numerical studies of the ISM with supernova explo- 
sions on a galaxy scale, we do not introduce simple 
evolutionary models for the SNR, such as the Sedov 
solution, and heating efficiency for the ISM due to 
a supernova. With our code, two-dimensional evolu- 
tion of blast waves caused by supernovae in an inho- 
mogeneous and turbulent medium with global rota- 
tion is simulated explicitly. Therefore, we can trace 
consistently the thermal and dynamical evolution of 
the ISM around the star forming regions and the as- 
sociated supernovae remnants and superbubbles (c.f. 
Norman fc Ikeuchi 1989| ) 

The hydrodynamic part of the basic equations 
is solved by the third-order AUSM (Advection Up- 
stream Splitting Method) flLiou fc Stcffcn 1993| ). Af- 
ter testing this code for various hydrodynamical 1-D 
and 2-D problems, we find that AUSM is as powerful 
a scheme for astrophysical p roblems as are the PPM 
( Woodward 



Colella 1984|) and Zeus ( [Btonc fc Nor 



i 1992| ) codes. More details about our numerical 



ma 

code and test results are described in Wada & Nor- 
man (2000, in preparation). 

We use 1024 2 Cartesian grid points covering a 8 
kpc x 8 kpc region. The spatial resolution is 7.8 
pc. A periodic Green function is used to calculate 
the self-gravity fo r the 8 kpc x 8 kpc region with 
2048 2 grid points (|Hockney fc Eastwood 198l|). The 



the cooling term in equation (3). 

The initial disk is an axisymmetric and rotation- 
ally supported (R = 3.7 kpc) with uniform surface 



density, T, g = 12M Q pc , an d the total gas mass is 
5 x 10 8 M Q (|Kim et al. 1998|). Random density and 



second-order leap-frog method is used for the time 
integration. We adopt implicit time integration for 



temperature fluctuations are added to the initial disk. 
Amplitude of the initial fluctuations is less than 5 % 
of the unperturbed values and have an approximately 
white noise distribution. The initial temperature is 
set to 10 4 K (R < 3.7 kpc) and 10 2 K (R > 3.7 kpc). 
The reason why we chose the low temperature in the 
outer region is to avoid the numerical artifact of the 
boundaries, i.e. reflection and generation of waves or 
shocks, for the initial evolution of the gas disk. In 
ghost zones at the boundaries of the calculating re- 
gion (i.e. 8x8 kpc 2 ), all physical quantities remain 
at their initial values during the calculations. From 
test runs we found that this boundary condition is 
much better than 'outflow' boundaries, because the 
latter cause strong unphysical reflection of waves at 
the boundaries. 

3. RESULTS 

3.1. Evolution and Structure of the ISM 

3.1.1. Model without Star Formation 

Figure 1 shows the time evolution of density and 
temperature of a model without star formation and its 
feedback (hereafter we call this model "model NSF"). 
Due to gravitational and thermal instability in the 
gas disk, clumpy fluctuations evolve in the first 10 8 
yr, and then the clumps merge and form larger struc- 
ture. They are deformed by the local tidal field and 
global shear in the non-linear phase, and as a re- 
sult, filamentary structure is formed. Higher density 
clumps (E 9 > 1O 3 M pc™ 2 ) are formed in the fila- 
ments due to the gravitational instability, or collisions 
between the filaments. The high temperature cavities 
(T > 10 5 K) are formed due to shock heating. The 
dynamical time scale for the large (L ~ 200 pc) high 
temperature (T ~ 10 5 K) region seen in the temper- 
ature panel at t = 800 Myr is Td yn ~ L/Av ~ 200 
pc/20 km s™ 1 ~ 10 7 yr. The cooling time, t coo i on 
the other hand, for this low density region (n ~ 10 -3 ) 
is ~ 2 x 10 7 yr (Spitzer 1977), which is comparable to 
Tdyn- In fact the high temperature cavity is a tentative 
structure that lasts ~ 10 7 yr. Photo-electric heating 
and UV irradiation contribute to form T <~ 10 4 K 
gas. In the high density (£ s > 1O 2 M pc™ 2 ) fila- 
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ments and clouds, the temperature is less than 100 K 
because the radiative cooling is effective. The global 
strucure of the disk does not change significantly after 
400 Myr, and it reaches a quasi-steady state. Figure 
2 is time evolution of the volume filling factor, /„, of 
each temperature level and time evolution of the max- 
imum density. This also represents that the global 
(> kpc scale) system reaches a quasi-stationally state 
after t ~ 300 Myr. This period is comparable to 
the local free fall time tff for the initial density, 
t ff ~ 2OO(S g /12M pc- 2 )- 1 / 2 Myr. At the global 
quasi-steady state, the warm gases (T g = 100 — 10 4 
K) occupy a large volume (~ 60 % of the whole calcu- 
lating region), and f v ~ 20 % for the cold (T g < 100 
K) and /„ - 20 % for hot gas (T g = 10 4 - 10 5 K). The 
maximum density in the same plot show slower evolu- 
tion after t ~ 300 Myr than the non-linear evolutional 
phase. 

One should note, however, that the globally stable 
state does not mean the local filamentary and porous 
structure whose scale is less than 1 kpc. It is rather 
quite dynamic and transient, which is similar to the 



past turbulent ISM models 


(Bania& Lyon 1980; Chi- 


anj & Prendergast 1985 


( 


jhiang & Bregman 1988]; 


Rosen & Bregman 1995 


; Vazquez-Semadeni ct al. 


19!! 5; Passot, Vazquez- Scmadcni, & Pouquct 1995|; 


Gazol & Passot 1999]). Clump mass is increasing dur- 



ing the initial linear and non-linear evolutional phase 
(t < 200 Myr) due to the mass accretion and merg- 
ing with other clumps or filaments, but in the quasi- 
stationally phase, disruption processes of the clumps, 
such as local tidal field (e.g. interaction between 
clumps, or shear) or shocks, prevent monotonic in- 
crease of the mass of each clump. Since we introduce 
the cutoff temperature for the cooling (i.e. 10 K), 
the gaseous pressure prevent the dense clouds further 
collapsing towards singularities. The angular momen- 
tum also supports the clouds (see §3.3). 

The density and temperature structures (Fig. 1) 
are similar to those in the model of Wada & Norman 
(1999), where a smaller disk is investigated (the ra- 
dius is 1 kpc), but the present model shows kpc-scale 
inhomogeneity and a more asymmetric distribution 
against the galactic center. This difference is prob- 
ably caused by the difference of the rotation curves 
used in the two models: a rigid rotation or differen- 
tial rotation. With a rigid rotation, random and tur- 
bulent motion dominates the circular rotation in the 
central region. With the global shear, on the other 
hand, global spirals develop toward the galactic cen- 



ter. 



3.1.2. Models with Star Formation 

Figure 3 represents density and temperature maps 
of the model SF2 at t = 833 Myr. The model SF1, 
which has ten times smaller threshold density for the 
star formation criterion (see §2), shows very similar 
density morphology and temperature structure. The 
red regions in the temperature map are hot gaseous 
regions where T g > 10 6 K. They are young (< 10 6 
yr) supernova remnants. Figure 3 also shows that 
most young supernova remnants are not axisymmet- 
ric. This is caused by that the background ISM is 
highly inhomogeous, and the radiative cooling at the 
dense filaments is so effective, and that it prevents 
the blastwaves expanding axi-symmetrically. Typical 
size of the hot cavities is less than 500 pc. Some bub- 
bles together form kpc scale 'super bubbles'. The hot 
region around (—3, —1) is one of the examples. Note 
that the size of the cavities is expected to change away 
from the galactic plane in a 3-D model (see also §4). 

Although the filamentary and clumpy structures of 
the gas which result in the star forming models are 
similar to those of model NSF, the kpc scale inhomo- 
geneity seen in the the model (Fig. 1) is not apparent 
in the SF models. The large scale inhomogeneity in 
the model SF1 and SF2 is less prominent than that 
of the model NSF. This is because that the time scale 
for making the kpc-scale holes, which is a dynamical 
time scale, ~ 10 s yr, is much longer than the the time 
scale for supernova explosions and evolution of the 
blastwaves, <; 10 6 yr. Approximately 10 3 supernovae 
per kpc 2 explode in this model during ~ 10 8 yr. The 
kpc-scale low density cavities, which is formed due 
to evacuation of gas by the effect of the gravitational 
and thermal instabilities, cannot evolve under such 
frequent supernova explosions. 

We find that the supernova rate is fluctuate in a 
time scale ~ 10 7 yr, but it stays in a range 4-10xl0 _4 
yr" 1 during 6xl0 8 yr in model SF2 (Fig. 4). This 
behavior also appears in model SF1, but the SN rate 
becomes smaller for the larger threshold density (~ 8- 
12xl0 -4 yr" 1 ). ROSAT observed 46 supernova rem- 



nants and more candidates in the LMC (Haberl & 
Pietsch 1999| ) . If we assume the ages of SNRs are 
between ~ 10000 yrs and ~ 30000 yrs, then we have 
SN explosion every 250 yrs or 750 yrs. This indi- 
cates about 0.0013/yr or 0.004/yr for SN rate in the 
LMC. If we use the larger threshold density for the 
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star formation in our model than that in model SF2, 
this results in smaller supernova rate than ~ 0.001 
yr _1 . Therefore the threshold density much larger 
than 400M© pc -2 would be excluded for modeling 
the LMC, if we assume that the star formation effi- 
ciency is 10 % or less and the standard IMF. 

In Figure 5, the star particles at t = 245 Myr 
in model SF1 are plotted, and one may compare it 
to the density distribution (the right panel) of the 
same snapshot. Massive stars are not uniformly dis- 
tributed, but they form clusters ('OB associations'). 
It is notable that the distribution of the star clusters 
does not necessarily correlate to the inhomogenous 
gaseous structure. In other words, cavities are not 
necessarily associate with the 'OB associations'. This 
is also seen in simulations by other groups, for exam- 
ple 2-D simulations of a highly compressible iso baric, 
no n-selfgravitatio nal fluid with star formation ( Scalo 
& ( phappcll 19991) . 



3.2. Line Emission Maps and Comparison with 
the Observations 

With the numerical simulations, we have density, 
temperature, and velocity fields (spatial resolution: 
7.8i pc). Using this information as an input, we com — 



putl o H I 21 cm brightness map and CO (J~l 0) lino 



map for model NSF and SF1 and SF2. These maps 
can be directly compared with the recent high reso- 
lution surveys of the LMC: H I with ATCA (Kim et 
al. 1998 ) and CO (J=l-0) with NANTEN QFukui ct 
il. JMQSL . 



Willi Ihc iimiieiical lesulls desuibed above, the 
following procedure is followed to compute a H I 21 
cm brightness map. It is assumed that atomic hydro- 
gen is in mostly neutral form in those regions where 
the temperature is at least a factor of 2 less than 8000 
K, a value typical of (mostly ionized) HII regions, and 
motivated by measurements of the kinetic tempera- 
tures of HI clouds and HI intercloud gas. It is also as- 
sumed that the H I level populations of interest follow 
a thermal distribution in those regions along the line 
of sight, and that this mostly neutral gas dominates 
the integrated emissivity. The line of sight is face- 
on, i.e., perpendicular to the grid used for the hydro- 
dynamic simulations. The two-dimensional density 
indeed is used as a column density in the radiative 
transfer calculations. The latter are done in three di- 
mensions by assuming a scale height for the neutral 
gas of H=100 pc, thus converting the column density 
in a (constant) local density for each point along the 



line of sight. It is this three-dimensional grid, with 
two-dimensional hydrodynamic information, that is 
used to determine the H I central brightness temper- 
ature in the Monte Carlo procedure (Spaans 1996), 
where the ambient radiation field and its interactions 
with matter is represented by a discrete number of 
"photon packages" . This method is three-dimensional 
and explicitly includes optical depth effects as well as 
the detailed velocity field of the hydrodynamical sim- 
ulations for photons that travel along lines of sight 
that are not face-on. These latter photon trajectories 
need to be incorporated when the level populations 
are not in thermal equilibrium, i.e., for CO. A local 
velocity dispersion AV of 1.29 x 10 4 T 1 / 2 cm s _1 is 
adopted for a kinetic temperature T, with a mini- 
mum of 0.5 km s _1 due to micro turbulent motions. 
For each grid point in the hydrodynamical simula- 
tion, the Monte Carlo radiative transfer then yields 
the integrated H I 21 cm intensity along the line of 
sight. 

For the CO (J=l-0) line, the approach is sim- 
ilar to the H I case, with the appropriate correc- 
tion factor in AV" for the different atomic weight of 
the CO molecule. Furthermore, a carbon chemistry 



is added to compute the abundance of CO (Spaans 



Van Dishoeck 1997), for a dust abundance equal 



to 1/5 of Solar. 



(van Dishoeck & Black 



This chemistry is well understood 
and care has been 



taken to include, for each computed line of sight 
and corresponding column density, the important self- 
shielding transitions of H2 and CO (c.f. Spaans et 



al. 1994). The ambient average interstellar radiation 



field, required for the ambient chemical balance, is 
determined by scaling with the LMC B band surface 
brightness in mag per square arcsecond with respect 
to Galactic. This typically yields an enhancement of 
a factor of 3 — 5 in the mean LMC energy density 
compared to Galactic. It is assumed, because the hy- 
drodynamical simulations are two-dimensional, that, 
to compute the local CO emissivity and the ambient 
chemical equilibrium, the gas density is given by the 
particular line of sight column density divided by a 
scale height of H = 100 pc for the neutral gas. The 
latter number H does not strongly influence the qual- 
itative features in the presented maps. 

Figure 6 presents H I 21 cm and CO (J=l-0) line 
brightness temperature maps calculated from the nu- 
merical data of the NSF model at t = 800 Myr, using 
the procedure described above. Figure 7 is the same 
plot as Fig. 6, but for model SF2 at t = 834 Myr. 
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Siz : of the H I filaments, shells, and holes in model Semadeni et al. 1995) ) . We find a rather weak but 



NSF is ~ 1 kpc. The filaments and cavities in the 
outer region (R > 2 kpc) of model SF2 are also kpc- 
scale. We would like to emphasize that large cavities 
(> kpc) and filaments in model NSF are NOT caused 
by a direct dynamical effect of star forming activities, 
but by the non-linear development of gravitational 
and thermal instabilities. 

The CO emission is localized in many clumps (size 
~ 10-100 pc) or clouds complexes (size ~ 0.1-1 kpc) 
in model NSF. Such CO "cloud" complexes could be 
sites for active star forming regions, such as the 30 
Dor star forming region in the LMC. Model SF2 shows 
more uniform distribution of "stars" than in model 
NSF. 

Figure 8 (a) and (b) are mass spectra of molecu- 
lar clouds obtained from the CO brightness temper- 
ature (Tb) distribution of models NSF and SF2. We 
identify clouds and their mass by the following pro- 
cedure. Assuming a threshold Tg for the CO map 
(Figs. 6 and 7), then we have many islands of CO 
emission. Most islands are not round in shape, but 
elongated or filament-like shape. We derive the mass 
of the islands using the surface density of the simula- 
tion data. Therefore the mass in Fig. 8 is the total 
gas mass of each island, not the virial mass. Here we 
plot three histograms for three different thresholds. 
The mass spectra shows a power-law, which is roughly 
dN c /dM c cx M" 1 - 7 for model NSF, where dN c is num- 
ber of clouds between the mass M c and M c + dM c , but 
this slope does not significantly depend on the thresh- 
old brightness. On the other hand, the spectrum of 
the star forming model steeper than the model with- 
out the energy feedback especially for small threshold: 
dN c /dM c cx M" 1 - 7 , M~ 2 ' 3 , and M" 2 ' 7 for model SF2 
with T B = 100, 50, and 30 K, respectively. The be- 
havior of the slope to the threshold Tg in the model 
NSF and SF2 means that the stellar energy feedback 
changes structure of low density envelope of the dense 
clumps. 

The NANTEN survey discovered about hundred 
molecular clouds in the LMC, and the mass spectrum 
of them is dN/dM vir cx M^ - 5±0A for a range between 



10 5_e M© QFukui et al. 1999| ), where M vir is the virial 
mass of the cloud estimated by using the observed line 
width. In our model, smaller clouds [M c < 10 5 AT Q ) 
tend to have smaller mass compared to their virial 
mass estimated from their internal velocity disper- 
sions. In other words, the smaller clouds are not in 



positive correlation between the cloud mass M c and 

1/2 

the virial mass Af vir in our model. The mass spec- 
trum of the model SF2 would be dN c /dM vir cx M~^ A , 
if we use the lowest threshold Tb — 30 K. A simi- 
lar discussion for the interpretation of the cloud mass 
spectra based on two-dimensional numerical simula- 
tions of the interstellar medium has been given by 
Vazquez-Semadeni et al. (1997) (see §4). 

The maximal mass of the clouds is approximately 
1O 6 ' S M and 10 5 - 5 ~ 6 Af Q in models NSF and models 
with the star formation (SF1 and SF2). This implies 
that very massive clouds (> 10 6 M Q ) are difficult to 
grow under the frequent supernova explosions. The 
maximal masses of the observed molecular clouds are 
about 1O 6 5 M QFukui et al. 1999|), which seems to 



prefer the model NSF. However, this does not mean 
that the star formation and its energy feedback are 
not important for shaping the ISM in the LMC, but 
it implies that the largest clouds could evolve in pref- 
erentially an environment where SNe are not frequent. 
Though one should be careful to make comparison be- 
tween the observed and numerical mass spectrum, be- 
cause the definition and identification of the "cloud" 
is not exactly identical, and it is related to noise level 
in observations. 

We also find that the H I map of the computa- 
tional model shows similar statistical properties to 
the observed H I map. Figure 9 shows the H I 
luminosity function, i.e. the histogram of the ob- 
served and numerical H I map of the LMC. The 
H I brightness temperature has been computed from 
Tb = S u \ 2 /2k B Q s b- S v is the Hi flux density, k B 
is the Boltzmann constant, and fl s b is the solid an- 
gle of the synthesized beam of ATCA mosaiced map. 
The H I intensity is determined from the integral of 
the brightness temperatures J Tsdv over the peak H I 
line profile, where dv is the channel width in kilome- 
ters per second. The H I luminosity functions appear 
to be log-normal like distribution for both observed 
and model NSF. The distributions of the H I lumi- 
nosity function of the model NSF and SF2 for lower 
surface brightness (Tb < 1000 K) arc similar to the 
observed one. The model SF2, on the other hand, 
shows excess above Tb > 1000 K than the observa- 
tions. These high brightness regions are originated 
from shock compressed, dense gas. 



equilibrium, but rather in a transient phase (Vazquez 
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3.3. Position- Velocity Diagrams 

Position- Velocity (PV) diagrams give us informa- 
tion on the kinematics of the ISM. In Fig. 10 (a), we 
show the PV diagram, in which y-componcnt of ve- 
locity (v y ) is integrated through x-positions, for the 
gas £ g < 1O 2 M pc~ 2 in model NSF. Many arc-like 
structures in Fig. 10 (a) look like expanding shells 
originating from explosional phenomena in the ISM. 
However, they are NOT caused by explosions because 
there is no energy input due to supernovae in model 
NSF. The arcs in the PV diagram are actually caused 
by the gases that form filament-like structure seen in 
the density map (Fig. 1). The filaments and shells 
exhibit non-circular motion of the order of 10 km s _1 , 
and their shape continuously changed. Here we would 
like to emphasize that, it is hard to distinguish, on 
the PV diagrams, between expanding shells and such 
shell-like structures which are changing in shape due 
to local random motion. 

From Fig. 10 (b), which is the same as Fig. 10(a), 
but for S s > 10 2 M Q pc~ 2 , we find that the dense 
and compact clumps, rotating with ~ 10 — 30 km s _1 , 
which are recognized as steep 'dotted lines' on the PV- 
diagram. Here we can identify about 80 clumps, and 
about half of them show retrograde rotation against 
the sense of galactic rotation. A representative one 
can be seen at (x, v y ) = (—0.5, —50) and less promi- 
nent one is at (1.5,25) and (3,75). Rotation of the 
clouds is important for the internal structure, motion, 
and star formation in the molecular clouds. Unfortu- 
nately, the spatial resolution in the NANTEN sur- 
vey (~ 40 pc) is insufficient to resolve the rotation 
of each molecular clouds of the LMC. High resolution 
observations and statistical analysis of kinematics of 
molecular clouds in external galaxies are necessary 
to understand the formation of molecular clouds and 
star formation processes. 

Figure 10 (c) shows the similar features shown in 
Fig. 10 (a), but for model SF1 at t = 610 Myr. Arc- 
like structures seen in the PV map for the NSF are 
not clearly shown in this map. This means that the 
linc-of-sight velocity field is much more random and 
chaotic than the NSF model, and there is no promi- 
nent large-scale coherent motion in the ISM of the 
model with stellar energy feedback. 

4. DISCUSSION 

The comparison of our model calculations with the 
observations of the LMC indicates that the star for- 



mation models are more consistent with the forma- 
tion of small (<C kpc ) and hot bubbles (T g > 10 7 
K), which have been detected by X-ray diffuse emis- 
sion. The power law slope of cloud mass spectrum of 
the models is close to that of observed one. On the 
other hand, the model calculations without the stel- 
lar feedback show even better agreement with global 
inhomogeneity of the HI and CO gas in the LMC, 
and also the HI distribution function is well repre- 
sented in the model NSF. Therefore, it is most likely 
that energy feedback to the ISM from massive stars 
and supernovae explosions are important processes for 
producing hot gas, but they do not necessarily dom- 
inate dynamics and formation processes for the large 
scale inhomogeneity observed in H I and CO gas in the 
LMC. However, one should note that the initial con- 
dition of the models (i.e. axisymmetric and uniform 
density distributions) causes nearly uniform star for- 
mations in the whole disk. If we begin the simulations 
from a more inhomogeous disk, stars will be formed 
non-uniformly, and effects of the stellar feedback are 
different, depending on location of the disk. We sus- 
pect that the more realistic model of the ISM in the 
LMC-type galaxy would be between our two extreme 
cases, but that can not be achieved by changing the 
threshold density of the star forming model. 

In our models with stellar feedback, typical sizes of 
the hot cavities are less than 500 pc. Some 'superbub- 
bles' can be formed in a outer disk region, and their 
linear sizes are ~ kpc (see Fig. 3). Rosen & Breg- 
man (1995) revealed in their two-dimensional, two- 
fluid (stars and gas) simulations in a disk galaxy that 
hot bubbles are formed by stellar activities, and their 
sizes depend on the energy injection rate. In their 
simulations, the linear size of the bubbles is up to 
- 500 - 1000 pc for the supernova rate is 0.0075-0.03 
yr _1 . This rate is 1-2 orders of magnitude larger than 
that of our star forming models. Local, three dimen- 
sional MHD simulations of the ISM with supernova 
explosions by Korpi et al. (1999) show that the lin- 
ear size of the superbubbles are typically 200-400 pc, 
which is consistent with our SF models. Gazol-Patino 
& Passot (1999) compute the evolution of the ISM 
in a region of the Galactic plane of size 1 kpc 2 in 
two-dimensional periodic domain. They found that 
superbubbles, and the largest one has linear scale is 
~ 500-1000 pc, due to about 3000 supernova explo- 
sions in ~ 10 7 yr. The background density of their 
simulation is equivalent to that in the outer region of 
our simulations. The local (< kpc scale) structure of 
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the ISM of our global simulations, where hot bubbles 
and cold filaments coexist, is also similar to the local, 
2-D simulations of the ISM with a periodic bound- 
ary condition (Vazquez-Semadeni et al. 1995). In a 
conclusion, the past local models with supernova ex- 
plosions, in which superbubble formation is reported, 
consistent with our global models with stellar energy 



feedback concerning the size ot the superbubblcs. 



ful to make direct comparison between results in the 
present paper and VBR97. 

Our result implies that there are two mechanisms 
of cavity formation. One is the pure evacuation 
of gas from the low-density regions by the effect of 
the gravitational and thermal instabilities. Similar 
evacuation phenomena (Elmcgrccn 1994; Vazquez- 



Semadeni, Passot, & Pouquct 1996) are also observed 



The cloud mass spectra of the model SF1 and SF2 
are steeper than that of the model NSF (Fig. 8). 
In other words, the stellar energy input affects the 
evolution of molecular clouds. Vazquez-Semadeni, 
Balleseros-Paredes, & Rodriguez (1997) (hereafter 
VBR97) analyzed their 2-D hydromagnetic simula- 
tions and showed that the mass spectra have the form 
dN c /dM c oc M~ 1,44±ai (c.f. their Fig. 3 and our Fig. 
8) , and they reported that the mass spectra from ob- 
servations are consistent with that from simulations 
in which the density filed had been shaped by stel- 
lar activity. Note that our mass spectrum of the star 
formation models, the slope is shallower for smaller 



clo jids [M c < W' 6 Mq) than that tor massive clouds 
(dfi c /dM oc M~ 2 ). Therefore our results might be 
consistent with the results of VBR97 (see also §3.2). 
However, direct comparison between these two results 
should not be straightforward, because there are num- 
ber of differences between our simulations and theirs, 
on the numerical scheme (AUSM vs. pseudo spec- 
tral method), the boundary condition (global simu- 
lations vs. local periodic boundary), and the cooling 
curve especially for T g < 100 K (A ^ vs. A = 0). 
The maximum density contrast is about 5000 in their 
model, but it is about 2 x 10 6 in our model. This 
difference is caused probably due to the difference of 
the cooling curve for the cold gas, and also due to 
the numerical scheme. In their numerical method, 
they added a mass diffusion term to the continuity 
equation in o rder to smooth out the density gradi- 
ents (see also Vazqucz-Scmadcni ct al. 1995 ). This 
could affect the structure of shocks and dense gas, i.e, 
the cloud mass spectrum. The recipe for the stellar 
energy feedback are also different. VBR97 assumed 
that once a star is formed, then it remains fixed with 
respect to the numerical grid. In our models, the star 
particles are orbited in the self-gravitational potential 
of the gas and the external fixed potential. The last 
point, however, would not be important, if number 
of stars is large enough, or the ISM is fully turbulent. 
VBR97, on the other hand, include the magnetic field 
in their models. Therefore, again, one should be care- 



in the simulations of the turbulent ISM without SNe 
(Passot, Vazquez-Semadeni, & Pouquct 1995). The 
other is the formation of cavities by the effect of SNe, 
and especially by the synchronized explosions of stars 
in OB associations. The two processes are essential 
to produce the whole structure of the ISM. The lat- 
ter is necessary for hot (T g > 10 6 K) component of 
the ISM. These two type of low-density regions are 
also observed in the two-dimensional, local simula- 
tions of the turbulent ISM with SNe by Gazol & Pas- 
sot (1999). 

As mentioned in §1, the origin of supergian t H I 



holes in galaxies has been controversial (see also Wal 



ter fc Brinks 1999| ) . If the observed kpc-scale holes 



and shells arc caused by cxplosional phenomena only, 
one needs highly energetic events, like Gamma Ray 
Bursts ( Loeb fc Perna 1995 ). However, our numeri- 
cal simulations give an alternative explanation about 
the origin of the large scale holes: the non-linear evo- 
lution of the multi-phase gas disk. In the multi-phase 
ISM, most of the gas mass is concentrated in the cold, 
dense clumps and filaments, but volume filling factor 
of such component is much smaller than that of hot, 
diffuse gas. The hot regions are surrounded by dense 
filaments as seen in Fig. 1 and Fig. 3. Therefore the 
multi-phase ISM in a quasi-steady state is naturally 
porous. Our results suggest that if the energy feed- 
back from massive stars are not effective, kpc-scale in- 
homogeneity can be evolved in a disk in about several 
10 8 yr. Dense filaments are changing in shape due to 
local random velocity field as well as the global shear, 
and often show kinematics features similar to those 
of "expanding shells" as seen in the position-velocity 
map (§3.3). We suspect that many supergiant holes 
and shells in dwarf galaxies do not have an explosion 
origin, but that supergiant shells far outside the disk 
plane can not be formed without cxplosional events. 

The work presented here yields the model for the 
ISM in a LMC-type galaxy. Nevertheless there are 
couple of things that one should consider for con- 
structing a complete numerical model of the LMC 
including the interaction with the Small Magellanic 
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Cloud, the optical bar, and non-uniform UV radia- 
tion field due to active star forming regions, such as 
the 30 Dor regions. The interactions with the SMC 
are important events for the star formation history in 



the LMC and the Magellanic Stream ( Murai fc Fuji 
mo,o 1980; Vallenari et al. 1996| ; pardiner fc Noguchi 
19£ 6|) . The off-center optical bar could also perturb 
the global structure of the ISM and star formation in 
the LMC ( pardiner, Turfus, fc Putman 1998j ). How- 
ever the H I mapping (Kim et al. 1998) does not show 
clear evidence of the stellar bar, namely off-set shocks 
which are often seen in barred galaxies. It would be 
interesting to investigate the effect of the off-center 
stellar bar in our model. This might contribute to 
the formation of the active star forming region, such 
as 30 Dor region. 

In the present paper, we have not solved the ver- 
tical structure of the ISM. It is expected that the 
hot component behaves differently in three dimen- 
sions. The hot gas, which is above 10 6 K, cannot be 
confined in t he disk plane ( Rosen fc Bregman 1995| ; 
Avillcz 1999| ), and it would probably been blown out 
from the disk plane. The scale-height of the hot gas 
should be larger than the cold gas, and, as a result, 
the volume filling factor of the hot gas would be larger 



aw av from the disk plane. Thus the radiative cooling 



would be less effective in the hot gas, in three di- 
mensions. This may affect interaction between cold 
and hot components, and the feedback process on the 
ISM. For example, it is more difficult to form the su- 
pergiant holes by SNe. In a subsequent paper, we 
will extend our method to three-dimensional model- 
ing, and investigate these problems. 

5. CONCLUSIONS 

Using high resolution hydrodynamical simulations, 
we have computed that the global dynamics and 
structure of the multi-phase ISM in a LMC-type 
galaxy. Due to gravitational and thermal instabil- 
ity in the gas disk, clumpy fluctuations evolve, and 
then the clumps merge and form filamentary struc- 
ture in the non-linear phase. Higher density clumps 
are formed in the filaments due to the gravitational 
instability, or collisions between the filaments. Our 
numerical model with the star formation allows us 
to provide the evolution of the blastwaves due to su- 
pcrnovae explosions in the rotating, inhomogeneous, 
multi-phase, and turbulent-like media. We find that 
the supernova rate in the model with stellar energy 



feedback is typically of the order of 0.001 yr _1 during 
several hundred Myrs, but fluctuates rapidly (time 
scale ~ 10 Myr) by a factor of three or four. The 
model also shows that kpc-scale low density cavities 



seen in the observed H I map (Kim et al. 1998) 



are difficult to be formed under frequent supernovae, 
but SNe are necessary to form hot bubbles where the 
gaseous temperature is greater than 10 6-7 K. Our 
result suggests that there are two possible causes of 
low density regions in the ISM. The kpc-scale inho- 
mogeneity and arcs can be formed as a natural con- 
sequence of non-linear evolution of the multi-phase 
interstellar medium in a LMC-type galaxy. We find 
in the PV-diagram of the numerical models that fil- 
amentary structure in the non-star forming model, 
which are caused by the gravitational instability, has 
similar kinematics to the structure formed by SN ex- 
plosions. The dense clouds are rotating, and about 
half of of them show retrograde rotation against the 
sense of galactic rotation. Using the Monte Carlo ra- 
diative transfer code, we have computed H I and CO 
brightness temperature distributions, and compared 
them with those from the recent observations. The 
CO cloud mass spectrum in the model with stellar en- 



ergy feedback is similar to the observed one (Fukui et 
al. 1999), but H I distribution function is well fitted 
by the model without SNe. Therefore we conclude 
that the small scale structure and dynamics of the 
ISM in the LMC is mainly affected by the stellar ac- 
tivities, but the gravitational instability significantly 
contributes to the global morphology and dynamics 
of the interstellar matter in a kpc-scale. 
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Fig. 1. — Evolution of density and temperature distri- 
butions of a model without star formation and energy 
feedback (model "NSF"). The density and tempera- 
ture are Log-scaled, and their units are Mq pc -2 and 
K, respectively. Time is shown at each panel in an 
unit, 10 s yr. Note: Although the gray scale is as- 
signed for 1000 — O.OIMq pc -2 in this plot, the maxi- 
mum density is much grater than 1000 M pc -2 (see 
Fig. 2). 

Fig. 2. — Evolution of the volume filling factor for 
four different phases, i.e. T < 100 K (filled circles), 
100 < T < 9000 K (open circles), 9000 < T < 10 5 K 
(open diamonds), and T > 10 5 K (filled diamonds). 
The maximum density of the gas is also plotted (thin 
solid line). 

Fig. 3. — Same as Fig. 1, but for a model with star 
formation and energy feedback (model "SF2" ) at t = 
834 Myr. 

Fig. 4. — Evolution of a total supernova rate in model 
SF2. 

Fig. 5. — Star particles and density distribution at 
t = 245 Myr in model SF1. 

Fig. 6. — Brightness temperature (Tb (K)) maps of 
H I (left) and CO (J=l-0) (right) computed from data 
of model NSF at t = 800 Myr. The intensity of H I is 
Log-scaled. 

Fig. 7. — Same as Fig. 6, but for model SF2 at t = 
834 Myr 

Fig. 8. — (a) Mass spectrum of CO clouds for model 
NSF for three thresholds of the brightness tempera- 
ture, Tb = 30 (solid line), 50 (dashed line), and 100 
K (dotted line). The thick line shows dN c /dM c cx 
M" 1 - 7 . (b) Same as (a), but for model SF2. 

Fig. 9. — Probability distribution function (pdf) of 
H I. Thick line shows pdf of model NSF (N(T B ) is 
the number of cells for the brightness temperature, 
Tb, and N is the total number of cells). The stars 
show the histogram from the H I synthesis observa- 
tion, with N 1 / 2 error bars, which is normalized for 
log(T B ) = 2.8 K of the model pdf. 

Fig. 10. — Position- Velocity maps for (a) low density 
gas in model NSF, (b) high density gas in model NSF, 
and (c) low density gas in model SF1. 
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